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1 .  INTRODUCTION 


In  general,  engineering  computations  serve  to  analyze 

Xsific  physical  phenomena  with  the  aim  of  reaching  decisions 
at  proposed  designs;  that  is,  for  instance,  the  expected 


behavlqr^of  a  construction,  its  safety  features,  and  similar 
concerns?^ For  a  general  assessment  of  the  structure,  quali¬ 
tative  results  about  the  displacements,  stress-distributions, 
etc.  are  needed  with  an  engineering  accuracy  of,  say  10-201. 
On  the  other  hand,  the  indicated  decisions  are  based  usually 
on  relatively  few  data  items  with  a  higher  accuracy,  say,  in 
the  range  of  2-5U,  such  as,  the  displacements  or  stresses  in 
certain  known  critical  points,  stress-intensity  factors,  etc., 
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In  order  to  achieve  these  diverse  objectives  the -desired 
>  computational  procedures  should  include  --  among  others-^  the 
following  features: 


(i)  Facilities  for  effective  and  reliable  estimation  of 
the  errors  of  the  computed  results  1  These  estimates 
should  correspond  to  the  specific  £ims  of  the  parti¬ 
cular  computation. 

^  N>  (d3o  Adaptively  controlled  computations^*  i^n  order  to 

achieve  the  desired  error  tolerances  with  a  minimal 
cost  and  to  reduce  the  preparatory  work  required  of 
the  user. 

(ijCj)  Post-processing  facilities,,  which  are  designed  to 

allow  for  the  extraction  of  high-accuracy  data  items 
of  the  type  needed  for  the  critical  decisions. 


For  the  solution  of  classes  of  linear  problems  by  means 
of  finite-element  approximations,  algorithms  for  all  three 
features  are  now  beginning  to  be  well  understood.  On  the 
other  hand,  the  nonlinear  case  is  as  yet  not  well  developed 
and  much  still  needs  to  be\done  before  satisfactory  algorithms 
for  all  three  features  can  be^onstructed. 

This  paper  is  not  intended  ^jpresent  an  overview  of  the 


work  in  this  area.  Instead  we  address  certain  specific 
questions  relating  to  the  theoretical  and  computational  aspects 

ujie  • 

of  the  three  features^  In  particular,  the  design  of  effective 

i. 

adaptive  procedures  requires  a  clear  understanding  of  the  ob¬ 
jectives  that  are  to  be  achieved. _ Accordingly,  in  Section  2 

we  discuss  some  new  results  aboul(  the  concepts  of  feedback  and 
adaptivity.  In  line  with  the  summary^ nature  of  this  paper. 


I 

1 


1 


3 


these  results  are  only  illustrated  for  a  one-dimensional , 
linear  model  problem.  Then  Section  3  addresses  their  ex¬ 
tension  to  the  more  general  two-dimensional,  linear,  elliptic 
class  of  problems  which  can  be  solved  by  the  adaptive  solver 
FEARS.  In  this  generality,  some  of  the  results  can  be  tested, 
as  yet,  only  experimentally.  In  Section  4  we  then  turn  to 
some  of  the  post-processing  approaches  for  linear  problems 
which  may  meet  the  demands  for  the  third  of  the  indicated  de¬ 
sirable  features. 

All  the  material  in  Sections  2  to  4  concerns  linear 
problems.  As  noted  before,  for  nonlinear  problems,  the 
corresponding  theory  is  by  far  not  so  well  developed.  In  part, 
this  is  due  to  the  fact  that  the  problem- formulations  differ 
considerably  between  the  two  cases.  In  particular,  in  practi¬ 
cal  nonlinear  problems  we  are  faced  usually  with  the  need  for 
considering  parametrized  systems  of  equations  for  which  the 
solutions  form  manifolds  in  the  space  of  the  state  and  para¬ 
meter  variables.  The  problem  then  becomes  the  computational 
determination  of  the  characteristic  features  of  such  solution 
manifolds.  This  is  the  topic  of  Section  5.  For  these  para¬ 
metrized  equations  it  is  not  always  readily  apparent  how  to 
define  the  discretization  errors  needed  for  the  first  of  our 
mentioned  features.  In  Section  6  we  present  some  recent  re¬ 
sults  about  a  priori  estimates  of  such  discretization  errors 
which  in  turn  lend  themselves  to  the  development  of  the  de- 

T. 

sired  a  posteriori  estimates.  Finally,  in  Section  7  we 

summarize  some  of  the  computational  approaches  which  combine 

* 

the  trace  of  paths  on  the  solution  manifolds  with  mesh-refine¬ 
ment  strategies  and  hence  which  may  provide  the  basis  for  the 


second  of  the  cited  facilities  in  the  nonlinear  case. 

\ 

2.  FEEDBACK  AND  ADAPTIVITY 

The  term  "adaptation",  or  any  of  its  cognates,  is  widely 
used  today  in  scientific  computing,  but  unfortunately  it  re¬ 
mains  a  rather  ill-defined  concept  which  is  often  used  very 
loosely.  Various  definitions  of  adaptive  processes  have  been 
given  in  the  literature  of  such  fields  as  automatic  control 
theory  (see  eg.  [22])  and  artificial  intelligence  (see  eg. 
[14]).  A  definition  which  applies  directly  to  many  of  the 
procedures  of  numerical  analysis  is  presented  in  [18].  With¬ 
out  entering  into  the  details  of  the  latter  definition,  we 
explain  in  this  section  the  relevant  concepts  of  feedback  and 
adaptivity  and  summarize  some  of  the  possible  theoretical 
results  for  the  case  of  a  simple  model  problem.  The  next 
section  then  addresses  the  extension  of  these  results  to  more 
general  problems. 

Consider  the  two-point  boundary  value  problem 


L[u]  =  -g^a(x)  +  b(x)u  - 

f ( x ) ,  x  e  I  ■  (0,1), 

(2.1) 

u(0)  *  u(l)  -  0 

(2.2) 

where 

0  <  a  <  a(x)  <_  a,  0  <_  b(x)  £ 

8,  X  £  i. 

(2.3) 

As  usual,  we  define  on  H1(I)  *  hP-cd  the  bilinear  form 

•B(u,v)  »  |  (au’v1  +  buv)dx  (2.4) 

I 

and  notu  that  the  norm  on  H1(I)  is  equivalent  with  the 


energy  norm  ||u]|E  -  (B(u,u))1/2.  The  weak  solution 
uQ  e  H1(I)  of  (2.1)/(2.2)  is  defined  by 

B(uq,v)  -  P(v)  =  J  fvdx,  V  v  s  H1(I).  (2.5) 

I 

Let  F  be  a  family  of  finite-dimensional  subspaces 
°1 

S  C  H  (I).  For  any  S  e  F  the  approximate  solution  u(S)  e  S 
of  (2.1)/(2.2)  is  specified  by 

B(u(S),v)  -  F(v),  ¥  v  e  S  (2.6) 

and  hence  the  optimal  error  achievable  with  subspaces  of  F 
of  a  fixed  dimension  d  >  0  is 

<Ku0,d)  -  inf {  |  |uQ-u(S)  |  |E;  S  e  F,  dim  S  -  d}.  (2.7) 

For  example,  we  obtain  a  set  F  ■  Fp  by  considering  all 
meshes 

A:  0  -*£<*£<  ...  <  *£(„  -  1  (2.8) 

and  specifying  S  e  Fp  as  the  space  of  all  C°-functions  on  I 
which  are  polynomials  of  degree  at  most  r  on  each  sub interval 
Ij  ■  CxJ_l5xJ],  i  ■  l,...,m(A),  of  A.  This  selection  leads 
to  the  classical  h-version  of  the  finite  element  method  for 
(2.1)/(2.2)  with  elements  of  order  r.  If  we  use  instead  a 
fixed  mesh  but  vary  the  order  of  the  polynomial  pieces 

separately  on  each  subinterval,  then  we  arrive  at  the  p-version 

\ 

of  the  method  (see  eg.  [8]). 

.A  feedback  system  for  (2.1)/(2.2)  on  some  given  set  F 
may  be  characterized  as  a  process  which  produces  sequences  (or 
trajectories)  of  subspaces  {S^}"  c  F.  This  is  accomplished 


by  some  transition  function  x  which  produces  for  any  index 
(time)  k  the  successor  sk+1  e  F  on  the  basis  of  infor¬ 
mation  about  the  sequence  up  to  k. 

For  a  precise  definition  of  this  concept  we  refer  to  [18]. 
Here  we  present  only  an  example  of  a  simple  transition 
function  x  ■  tp  for  the  case  F  *  Fp.  For  any  S  ■  S(A)  t  Fp 
consider  for  each  1,  1  <  i  <_  m(A),  the  exact  (weak) 
solutions  z 1  of  (2.1)  on  1^  ■  Lx^^jX^]*  subject  to  the 
boundary  conditions  *  u(S)(x^_1),  z^x^)  *  u(S)(x^). 

Then  the  error  indicators 


ij  -  I |z1-u(S(A)) | lE>IA,  i  -  l,...,m(A), 


(2.9) 


can  be  computed  approximately  or  estimated  from  above  and 


below  in  terms  of  the  residuals  r£  »  L[u(S)]  -  f  on  I 
Moreover 


i* 


m(A) 


«(*)  -  t  T  (nj)2]1/2 

i-1  X 


(2.10) 


represents  an  estimator  for  the  error  | |uQ-u(S) | | p,  (see  [4], 
[5],  [6]).  Now  a  simple  mesh-refinement  strategy  constructs 
from  A  a  new  mesh  A'  by  halving  the  subinterval  1^  for 
which  is  maximal.  This  corresponds  to  the  transition 

function  x  ■  xp  which  produces  as  the  successor  x(S(A))  of 
the  current  space  S(A)  e  Fp  the  space  S’  •  S(A’)  e  Fp. 

,Now  any  feedback  system  on  a  set  F  will  be  called  an 
adaptive  system  if  it  is  optimal  with  respect  to  some  per¬ 
formance  measure,  or,  more  precisely,  if  under  this  measure 
its  performance  is  not  worse  than  that  of  any  other  feedback 
system. 


Generally,  a  performance  measure  is  a  function  u  which 


associates  with  any  trajectory  (Sk>  C.  F  produced  by  a  feed¬ 
back  system  an  element  y{Sk>  of  a  given  partially  order  set 
M.  In  the  simplest  case  we  might  use  the  set  M  ■  {0,1} 
which  allows  us  to  distinguish  only  between  good  (*0)  or  bad 
(«1)  performance.  Two  typical  performance  measures  of  this 
type  may  be  defined  as  follows: 


/" 0  if  lim  | |uQ-u(Sk) I lE  -  0 


u{sk}  -  { 


(2.11) 


V  1  otherwise 


u(S.  } 


*"l. 


0  if  | |uQ-u(Sk) | |£  <  C  ♦(u0,dim(Sk)), V  k  >  0 

(2.12) 


otherwise 


where  in  (2.12)  the  constant  C  depends  on  the  problem  and 

the  solution  u  -  but  not  on  k. 
o  ’ 

These  and  other  measures  are  discussed  in  [l],  [2],  [9]. 
We  mention  here  only  a  few  partial  results  for  our  model 
problem  phrased  not  necessarily  in  their  most  general  form: 


°1 

Theorem  2.1:  For  the  problem  (2.1)— (2.3)  and  any  uQ  e  H  (I) 

consider  the  simple  feedback  system  on  F  ■  Fp  defined  by  the 

transition  function  xp.  Then  the  resulting  trajectories 

{S.  >  satisfy  lim  | |u  -u(S.  ) | |_  ■  0  and  hence  the  system  is 
K  k-*»  °  k  a 

adaptive  with  respect  to  the  performance  measure  (2.11). 

For  the  performance  measure  (2.12)  more  restrictive 
conditions  are  required.  In  essence,  we  compare  here  the 


convergence  rate  for  the  computed  sequence  of  meshes  with  an 

"ideal"  rate.  For  this  we  need  to  know  —  irrespective  of  any 

feedback  system  —  the  existence  of  a  sequence  of  meshes  for 

which  some  comparison  with  the  ideal  rate  applies.  We  call 

this  a  comparison  sequence;  its  existence  is  a  condition  on 

°1 

the  solution  uQ  e  H  (I)  which  we  wish  to  approximate.  In 
essence,  a  theorem  of  the  following  form  can  then  be  proved 
(see  L9]). 

Theorem  2.2:  Suppose  that  the  coefficients  a,b  of  (2.1/3)  are 

°1 

sufficiently  smooth  and  that  for  the  given  uQ  e  H  (I)  there 
exists  a  comparison  sequence  of  meshes  with  certain  properties. 
Then,  for  the  feedback  system  on  Fp  defined  by  tp,  the  re¬ 
sulting  trajectories  satisfy 

I luo~u*Sk) I lE  I  c  <fr(uo,dim  S(uk)),  k  >  0 

with  a  constant  that  depends  on  the  problem  and  the  comparison 
sequence,  but  not  on  k.  Hence,  the  feedback  system  is 
adaptive  under  the  performance  measure  (2.12). 

The  result  implies  that  under  the  particular  conditions 
the  rate  of  convergence  is  independent  of  the  smoothness  of 
the  solution.  The  conditions  of  the  theorem  are  not  overly 

a  °1 

stringent.  For  example,  they  hold  for  uQ  •  x  -  x  e  H'(l), 
o  >  j  ,  in  which  case  the  rate  of  convergence  is  of  the  order 
of  (dim  S(Ak))“r  independent  of  a,  (see  [9]).  We  refer  also 
to  [5]  for  related  results. 

The  simple  transition  function  Tp  used  in  both  theorems 
is  computationally  co  tly  since  it  modifies  the  meshes  very 
little  from  t  step.  Moreover,  in  practice  we  can  work 


only  with  a  finite  segment  S1,...,Sn  of  any  trajectory.  In 
that  case,  it  is  reasonable  to  demand  that  the  computation  of 
the  final  approximate  solution  u(Sn)  is  at  most,  say,  twice 
or  thrice  the  cummulative  cost  of  the  computation  up  to  that 
point . 

For  this  we  have  to  introduce  an  appropriate  stopping 
criterion  which,  obviously,  needs  to  be  based  on  an  error 
estimator  such  as  (2.10).  The  construction  of  such  a  stopping 
rule  is  suggested  by  the  following  two  results: 

°1 

Theorem  2.3:  For  uQ  e  H  (I)  and  under  the  condition  (2.3) 
we  have 

C2  e ( A)  <  | |uq-u(S(A) ) | |e  <  Cx  e(A) 
where  the  constants  ci»C2  >  0  are  lndePendent  of  S(A). 

Theorem  2.4:  Under  the  conditions  of  theorem  2.2  we  have 


For  some  results  about  the  adaptivity  of  feedback  systems 
with  transition  functions  incorporating  such  stopping  criteria 
see  [l],  [9].  It  should  be  noted  also  that  the  results  in¬ 
dicated  here  remain  valid  for  finite  sets  of  right-hand  sides 
f  in  (2.1). 

3.  MORE  GENERAL  PROBLEMS  AND  THE  FEARS  SYSTEM 

The  results  sketched  in  the  previous  section  allow 


10 


extensions  to  more  general  problems  in  two  space  dimensions. 

In  particular,  results  presented  in  [2]  address  the  family  of 
problems  and  feedback  approaches  upon  which  the  FEARS-program 
[Finite  Element  Adaptive  Research  Solver]  is  based. 

For  a  detailed  description  of  FEARS  we  refer  to  [16J  and 
[23].  In  brief,  FEARS  solves  certain  systems  of  elliptic 
equations  in  two  space  dimensions  and  permits  also  combi¬ 
nations  with  one-dimensional  problems,  as,  for  instance,  those 
needed  for  modelling  stiffeners  in  elasticity  problems.  Curved 
quadrilateral  elements  of  first  order  are  used.  As  in  (2.9) 
error  Indicators  are  associated  with  each  element  and  from 
them  an  estimator  (2.10)  of  the  solution-error  is  obtained. 

For  details  about  these  a  posteriori  estimates  see  [2], 

FEARS  utilizes  various  modes  of  mesh-refinement.  In 
particular,  a  so-called  long  pass  involves  the  recomputation 
of  the  solution  on  the  new  mesh,  while  in  a  short  pass  no 
such  recomputation  Is  performed  but  instead  decisions  about 
further  refinements  are  based  on  an  extrapolation  of  the  re¬ 
finement  patterns  observed  In  previous  meshes.  Thus,  by 
combining  these  short  and  long  passes  in  various  ways  differ¬ 
ent  feedback  mechanisms  are  obtained. 

Theorem  2.1  carries  over,  in  full  generality,  to  tra¬ 
jectories  generated  by  long  passes  In  FEARS,  and  the  same  is 
true  for  Theorem  2.3  (see  [3],  [**]).  On  the  other  hand,  the 
analogue  of  Theorem  2.2  is  not  available  in  sufficient 
generality.  Yet,  all  numerical  experiments  performed  with 
FEARS  suggest  that  such  a  theorem  should  be  valid  also  in 
this  case.  Theorem  2.H  does  extend,  however,  under  certain 
additional  assumptions  about  the  meshes  (see  [2]).  By 


11 


appropriate  use  of  the  short  passes  it  is  possible  to  ensure 
that  the  cost  of  the  computation  of  the  solution  on  the  final 
meshes  of  the  refinement  sequence  is  proportional  to  the  cumu¬ 
lative  cost  of  the  computation  up  to  that  point.  In  fact,  ex¬ 
perience  shows  that  the  proportionality  factor  is  of  the  order 
of  1.5  to  3.0. 

As  an  illustration  we  present  here  an  example  of  the  typi¬ 
cal  performance  of  FEARS.  More  specifically  we  consider  the 
problem  Au  *  0  on  the  domain  shown  in  Figure  1  with  the  indi¬ 
cated  boundary  conditions.  The  figure  shows  also  the  dotted 

8  *2 

Figure  1 


partition  of  the  domain  used  in  the  mesh-construction  of  FEARS 
(see  [23]).  The  solution  has  a  strong  singularity  at  the  ori- 

1/4 

gin  of  the  type  r  ;  that  is,  the  solution  belongs  to 

for  any  e  >  0.  Hence  for  sequences  of  uniform 
meshes  the  rate  of  convergence  is  of  the  order  of  N“  '  with 
respect  to  the  degrees  of  freedom  N.  This  is  certainly  a  slow 
rate. 

Figure  2  depicts  the  achieved  accuracy  —  measured  in  the 
energy  norm  —  as  function  of  the  number  of  elements  N*  (—  N). 
More  specifically,  the  true  errors  are  given  for  a  sequence 
(a)  of  uniform  meshes,  (b)  of  adaptively  constructed  meshes 


based  on  the  combination  of  long  (L)  and  short  (S)  passes 
LSLSSLSSSLSSSSLSSSL,  and  (c)  of  adaptively  constructed  meshes 
using  long  passes  only.  In  addition,  the  computed  error  esti¬ 
mates  are  shown.  All  data  are  in  percent  of  the  square  root 
of  the  exact  energy 


°  uniform  mesh 
>  adaptive  long 
passes  only 

o  adaptive 
combination  of 
long  and  short 
passes 

♦  error  estimator 


Figure  2 


The  figure  indicates  that  the  procedure  is  convergent  in  ac¬ 
cordance  with  the  two-dimensional  analogue  of  Theorem  2.1.  The 

rate  of  convergence  for  the  adaptively  constructed  meshes  has 

-1  /2 

the  maximum  rate  N  '  ,  which  confirms  experimentally  the 

analogue  of  Theorem  2.2.  Actually  we  see  a  slightly  better 
rate  because  the  higher  errors  are  here  principally  confined  to 
a  neighborhood  of  the  origin.  On  the  other  hand,  the  rate  of 
convergence  for  the  uniform  meshes  is  very  close  to  N"  '  .  In 
fact,  the  accuracy  of  8.5%  obtained  by  the  adaptive  mesh  of  537 
elements  would  require  a  uniform  mesh  of  more  than  10^  elements. 

The  adaptive  meshes  based  on  long  passes  give  better  re¬ 
sults  than  those  involving  a  combination  of  passes.  In  de¬ 
pendence  of  N  this  is  certainly  expected.  The  situation  is 
different  in  Figure  3  where  the  errors  are  shown  in  relation 
to  the  total  machine  time.  Here  the  effectiveness  of  the 
combination  of  long  and  short  passes  shows  up  very  clearly. 


TOTAL  MACHINE  TIME 


Figure  4  gives  the  earlier  mentioned  factor  of  the  cost  of  the 
final  long  pass  to  the  total  time  for  the  sequence  of  long  and 
short  passes.  Evidently  the  factor  is  of  the  order  of  two. 


RELATIVE  ERROR  IN  ENERGY  NORM(%) 


Figure  ** 


Finally,  Figure  5  shows  the  effectivity  index 
0  ■  est.  error/true  error  as  a  function  of  the  achieved  accu¬ 
racy.  Experience  shows  that  often,  but  not  always,  the  esti¬ 
mated  errors  are  smaller  than  the  true  errors.  But  for  the 
accuracies  of  10X  or  better  the  effectivity  index  exceeds  0.9. 
In  line  with  this,  experience  has  shown  that  the  stopping 
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RELATIVE  ERROR  IN  ENERGY  NORM  (%) 


criterion  based  on  the  estimators  is  very  reliable.  Generally 
the  effectivity  index  is  better  for  smoother  solutions;  it  is 
also  better  for  the  procedure  based  only  on  long  passes. 

The  results  discussed  here  for  this  single  example  are 
consistent  with  and  typical  of  all  the  experimental  evidence 
we  have  obtained  so  far  with  PEARS. 


4.  POST-PROCESSING 

As  mentioned  in  the  Introduction,  in  practical  engineer¬ 
ing  computations  Important  decisions  are  usually  based  on 
relatively  few,  but  highly  accurate  data  items.  A  very  ef¬ 
fective  approach  to  the  computation  of  such  data  is  the  appli¬ 
cation  of  post-processing  techniques  to  the  computer  so¬ 
lutions.  For  a  detailed  study  of  this  approach  we  refer  to 
[3].  Here  we  sketch  only  the  principal  ideas  for  the  simple 


model  problem 


-Au  *1,  on  n  ■  (0,1)  *  (0,1) 


(4.1) 


u  *  o,  on  an. 


(4.2) 


Suppose  that  we  are  interested  in  the  values  ^  «  u(0,0)  and 
♦2  *  ^  (1*0).  The  main  idea  is  to  write 

*i  -  !  u^dx  +  R±.  i  »  1,2  (4.3) 

n 

where  uQ  again  denotes  the  exact  solution  of  (4.l)/(4.2)  and 
the  Ri  are  integrals  that  may  be  computed  from  the  input 
data.  For  instance,  for  ^  we  may  use 

■  -Aip,  Ri  “  }  ^dx 

n 

ip(x)  ■  ~  (lu  |  j  x  |  1 2  -  luCjU+x^Hl+x2)]3^2). 

Once  (4.3)  is  available,  an  approximate  value  of  ^  can  be 
computed  by  replacing  in  (4.3)  the  exact  solution  uQ  with 
the  computed  solution  u(S)  (see  the  notation  of  Section  2). 
Then  the  error  can  be  written  in  the  form 


♦  l  -  #i  ■  B(u(S)-uQ,ip1(S)-ipol) 

where 


B( 


„\  _  f  r  3u  3v  .  3u  3v 

u>v)  ■  j  C3^  55^  *  jxj 


and  ipQl  is  the  solution  of  -Aip  *  subject  to  the 
boundary  conditions  (4.2)  while  ^^(S)  is  the  corresponding 


approximate  solution.  The  computation  of  ^(S)  can  be  done 
simultaneously  with  the  computation  of  u(S)  at  no  great 
additional  expense. 

Table  2  shows  a  comparison  of  the  results  for  4> ^  and 
♦2  obtained  both  from  the  computed  solution  u(S)  and  the 
post-processing  approach. 


No.  of  elements 
In  quarter  domain 


16 


6 


_ 

Effectlvity  index 

for 

1.07 

.98 

• 

VO 

VO 

Effectlvity  index 
for  $2 

.93 

1.02 

1.06 

Table  3 

5.  PROBLEM  FORMULATIONS  IN  THE  NONLINEAR  CASE 

By  far  the  majority  of  todays  finite -element  computations 
Involve  linear  problems;  but.  In  recent  years,  the  nonlinear 
analysis  of  structural  systems  has  become  Increasingly  im¬ 
portant  In  engineering  applications.  Not  unexpectedly,  there 
are  many  fundamental  differences  between  linear  and  nonlinear 
problems,  and.  In  fact,  the  problem  formulations  for  the 
latter  show  many  special  features  not  present  In  the  linear 
case. 

First  and  foremost.  It  has  to  be  recognized  that  the 
theoretical  approaches  and  computationally  techniques  depend 
strongly  on  the  type  and  properties  of  the  nonlinearities.  In 
structural  problems  we  may  distinguish  four  principal  sources 
of  nonllnearltles,  namely  (see  [11]) 

(1)  geometric  nonllnearltles.  due  to  nonlinear  strain- 
displacement  relations 

(11)  material  nonllnearltles.  due  to  nonlinear  consti¬ 
tutive  equations 


(iii)  force  nonllnearltles .  due  to  nonlinear  stress 
boundary  conditions 

(iv)  kinematic  constraint  nonlinearities,  due  to  non¬ 
linear  displacement  boundary  conditions. 

The  source  of  nonlinearity  affects  the  form  of  the  resulting 
nonlinear  equations  and  hence  the  solution  approaches.  The 
corresponding  numerical  problems  are  as  yet  not  completely 
understood.  This  is  particularly  true  for  the  nonlinearities 
of  type  (ii)  to  (iv). 

As  indicated  in  the  Introduction,  engineering  compu¬ 
tations  serve  to  predict  the  behavior  of  a  physical  system 
under  study  so  as  to  reach  decisions  about  proposed  designs. 

In  the  case  of  structural  problems,  the  term  "behavior"  refers, 
say,  to  deformations  in  response  to  the  action  of  external  in¬ 
fluence  quantities,  such  as  loads,  changes  in  material  pro¬ 
perties  and  geometrical  data,  etc.  This  means  that  the 
equations  of  the  system  also  depend  on  a  set  of  influence 
parameters  and  the  objective  Is  to  assess  the  behavior  of  the 
solutions  under  variations  of  the  parameters.  In  the  linear 
case,  this  variation  is  also  assumed  to  be  linear  and  we  need 
only  compute  a  few  specific  solutions  to  achieve  our  objective. 
In  nonlinear  problems,  however,  the  computation  of  a  few 
specific  solutions  provides  little  or  no  insight  into  the 
behavior  of  the  system  and  we  are  led  to  consider  the  set  of 
all  solutions  depending  on  all  the  parameters  in  a  specific 
range. 

Mathematically,  this  means  that  we  are  faced  with  para¬ 
metrized  equations  of  the  form 
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F(y,A)  -  b  (5.1) 

where  y  varies  in  a  space  Y  of  state  (behavior)  variables, 

X  in  a  parameter  space  A,  P  is  a  given  mapping  from  Y  x  a 
into  Y,  and  b  e  Y  a  prescribed  element  in  the  range  of  F. 
In  general,  the  set  of  all  solutions  (y,A)  e  Y  *  A  of  (5.1) 
forms  a  manifold  in  Y  *  A.  In  structural  mechanics,  this  set 
is  now  usually  called  the  equilibrium  manifold  of  the 
structure.  The  objective  then  becomes  a  computational  analysis 
of  the  form  and  characteristic  features  of  this  manifold. 

Clearly,  the  specific  features  to  be  determined  for  a 
particular  problem  depend  on  the  objectives  of  the  compu¬ 
tation.  Generally,  interest  centers  on  computing  solutions  on 
the  manifold  along  specified  paths  which  represent,  for  in¬ 
stance,  particular  loading-regimes,  etc.  Then,  the  boundaries 
of  the  stability  regions  on  the  manifold  are  to  be  explored, 
and  the  form  of  the  manifold  near  specific  singularities  is  to 
be  determined. 

Besides  these  general  objectives,  there  are  various 
specific  solution  data  that  are  needed  for  certain  decisions 
about  the  structure  under  consideration.  A  typical  example  of 
this  is  the  so-called  question  about  conservative  input  data. 
In  general,  the  external  input  data  influencing  a  structure 
are  only  known  within  certain  ranges;  for  Instance,  this  is 
certainly  true  for  the  loads  to  which  the  structure  will  be 
subjected  during  its  lifetime.  Thus,  for  the  design  it  is 
essential  to  know  the  values  of  the  parameters  which  represent 
worst  possible  cases.  These  are  called  conservative  input 
data.  Their  determination  represents  a  search  problem  on  the 
solution  manifold. 


2 


As  an  example,  consider  a  clamped,  thin,  shallow  circular 
arch  which  has  been  used  as  a  test  case  by  various  authors 
(see  eg.  [15],  [21]).  Let  U  and  W  be  the  radial  and  axial 
displacements,  R  the  arch-radius,  A  Its  cross-sectional 
area,  h  the  thickness,  and  E  Young's  modulus.  With  the 
dimensionless  displacements  u  ■  U/h,  w  -  W/h,  the  total 
potential  energy  —  non-dimensionalized  by  dividing  by 
EAR(h/R)2  —  Is  given  by 


\  1 0[(3?-u)  ♦  H  <a?)2}2'i9  ♦  tt  1  °(?rt)2<16  -  °2 1  °pua6- 


Here  p  -  p(6)  Is  the  dimensionless  radial  load  and  a^,a2 
are  dimensionless  constants.  Each  end  Is  assumed  to  be 
pinned;  that  Is,  we  have  the  boundary  conditions 


u<+8  )  -  0,  w( +0  )  -  0,  (+8  )  -  0. 

-  °  -  °  de2  ■  ° 


The  finite  element  approximations  Introduced  in  [25]  were 
applied.  More  specifically,  we  used  a  uniform  mesh  with  eight 
elements,  0Q  ■  15°,  and  the  constants  ■  3.8716  x  10“**, 

<*2  »  8.2752  x  io  corresponding  to  data  In  [15]. 

Suppose  that  the  following  two-parameter  family  of  loads 
Is  to  be  considered 


p(y,v)(9) 


(1  -  ^-(v-6)  )y,  for  max(-60,v— jp)  <  8  <  v 
o 


(1  -  ^-(0-v)  )y ,  for  v  <  0  min(0o,v+-jp) 
o 


0 


9 


otherwise 


Thus  the  load  Is  a  plecewlse-linear  hat-function  which  has  the 
value  v  at  6  »  v  and  is  zero  outside  the  interval  centered 
at  v  of  width  0q/2.  For  any  fixed  position  v  the  path  on 
the  manifold  determined  by  varying  y  encounters  a  first 
limit  point  which  represents  a  point  where  the  arch  is  expect¬ 
ed  to  buckle.  Figure  6  shows  the  variation  of  this  buckling 
load  with  the  value  of  v.  We  see  that  the  most  dangerous 
loads  occur  about  for  v  ■  +O.l60o;  hence,  these  values  repre¬ 
sent  the  desired  conservative  data.  The  computations  were 
performed  with  the  limit-point  facilities  of  the  continuation 
package  PITCON  (see  [19]). 


Figure  6 


6.  ERROR  ESTIMATES 


As  in  the  linear  case,  a  desirable  feature  for  nonlinear 
structural  computations  is  the  estimation  of  the  errors  of  the 
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computed  results.  But,  not  unexpectedly,  for  nonlinear 
problems  the  situation  is  considerably  more  difficult  and  much 
remains  to  be  done.  In  fact,  for  the  solutions  of  the  para¬ 
metrized  equations  of  the  previous  section,  it  is  not  always 

readily  apparent  how  to  define  the  error.  In  fact,  if,  say, 

$ 

M  denotes  the  equilibrium  manifold  of  a  given  Infinite- 
dimensional  equation  (5.1)  and  M  is  the  solution  manifold  of 
a  finite-dimensional  approximation  of  this  equation,  then  it  is 

me 

not  obvious  how  to  relate  points  of  M  and  M .  Frequently, 
the  parameter  space  A  of  (5.1)  is  finite-dimensional  and 
hence  there  is  no  need  to  discretize  X  e  A.  In  that  case,  we 
might  wish  to  compare  points  of  M  and  M  corresponding  to 
the  same  values  of  A.  But  already  simple  examples  show  that 
for  specific  A  there  may  exist  a  point  only  on  one  of  the  two 
manifolds. 

In  order  to  resolve  this  problem  it  is  necessary  to  con¬ 
sider  first  the  equation  of  a  priori  estimates  for  parametrized 
equations.  For  mildly  nonlinear  elliptic  boundary  value 
problems  and  one-dimensional  parameter  spaces  such  estimates 
were  first  developed  in  [10] .  A  different  approach  was  used  in 
[12]  to  provide  similar  results  for  more  general  boundary  value 
problems.  In  [13]  this  approach  was  extended  to  the  case  of 
any  finite-dimensional  parameter  space.  We  present  here  only  a 
brief  summary  of  the  results  in  [12],  [13]. 

Let  X,Y  be  Banach-spaces  and  F:  S  e  x  ♦  Y  a  mapping  of 
class  Cr(S),  r  ^  1,  on  an  cpen  connected  subset  S  of  X. 
Suppose  that  F  is  a  Fredholm  operator  of  index  m  :>  1  on  S 
and  that  the  regularity  set  R(F)  ■  (x  e  S;  rge  DF(x)  ■  Y)  is 
non-empty.  Then  R(F)  is  an  open  subset  of  X  and  for  any 


regular  value  b  e  F(R(F))  the  regular  solution  set 
M(b)  -  {x  e  R(F) ;  F(x)  -  b} 

in  an  m-dlmenslonal  Cr-manlfold  in  X.  For  proofs  we  refer  to 

[12]. 

For  any  theoretical  and  computational  study  of  M(b)  we 
need  to  specify  which  local  parametrlzatlons  of  the  manifold 
we  will  use.  At  a  point  xQ  e  M(b)  a  local  parametrlzatlon 


of  M(b) 

shall  be  a  triple  {V,A,T> 

consisting  of 

(1) 

a  closed  subspace  V  of 

X  such 

that 

VD  ker  DF(x_)  -  {0}, 

(11) 

an  isomorphism  A  e  L(Y,V) 

from 

Y  onto  V, 

(ill) 

an  m-dlmenslonal  subspace. 

T  of 

X  such  that 

X  -  V  ©  T. 


The  meaning  of  this  concept  Is  contained  In  the  following  re¬ 
sult  (see  [12]) : 

Theorem  6.1:  Under  the  stated  conditions  for  F  and  M(b), 
let  {V,A,T}  be  a  local  parametrlzatlon  of  the  manifold  at 
xQ  e  M(b).  Then  there  exists  an  open  ball  Be  T  with  0  e  B 
an  open  neighborhood  U  C  X  of  xQ,  and  a  unique  Cr-function 
n:  B  Y  such  that 

M(b)A  U  •  (x  e  X;  x  ■  x  +  t  +  An(t),  V  t  e  B}. 

o 

The  geometric  Interpretation  is  sketched  In  the  composite 
Figure  7.  In  many  applications  certain  parameters  are  identi¬ 
fied  explicitly;  that  Is,  one  has  a  natural  splitting 
X  ■  W  ©  Z,  dim  Z  ■  m,  where  W  Is  Isomorphic  with  the  range 


space  Y.  The  question  then  arises  at  what  points  xQ  e  W(b) 
this  natural  splitting  gives  rise  to  a  local  pararaetrization; 
that  is,  when  the  given  parameters  can  be  used  to  parametrize 
the  manifold  near  xQ.  Clearly,  the  condition  here  is 
W  n  ker  DF(xq)  ■  {0}  and  this  holds  exactly  if  rge  D^FCXq)  »Y. 
Then  {W,DwF(xo)“1,Z}  indeed  is  a  local  parametrization  of 
M(b)  at  xQ.  Such  points  are  called  nonsingular  points  of 
M(b)  with  respect  to  the  natural  parameters. 

At  points  xQ  e  M  (b)  where  WO  ker  DF ( xQ )  W  {0}  the 
natural  parameters  cannot  be  used  to  parametrize  M(b).  Such 
points  can  be  characterized  by  the  dimension  of  this  inter¬ 
section.  The  simplest  case  arises  when  dim(W  r\  ker  DF(xq))*1; 
these  are  limit  points  with  respect  to  the  direction  of  the 
intersection.  For  more  details  we  refer  to  [13]. 

We  turn  now  to  a  discretization  of  our  problem.  Once 
again,  it  is  useful  to  assume  the  existence  of  a  natural 
splitting  X  ■  W  ©  Z,  dim  Z  ■  m,  where  W  is  isomorphic  with 
Y.  More  specifically,  we  assume  that  there  is  an  operator  Q 
such  that 

Q  e  L(X,Y),  Z  -  ker  Q,  rge  Q|W  -  Y. 

Often  we  have  F:  S  CL  Y  x  Rm  -*•  y  In  which  case  Q  can  be 
taken  as  the  natural  projection  from  Y  x  Rm  onto  Y. 

The  natural  splitting  has  to  be  considered  since  in  such 
applications  only  the  state  space  W  is  discretized  and  the 
parameters  are  left  untouched.  Suppose  now  that  the  discreti¬ 
zation  is  defined  by  a  family  of  linear  maps  Ph  e  L(Y)  of 
finite  rank  indexed  by  a  positive  h  >  0  such  that 

lim  PKy  -  y  for  all  y  e  Y.  Then  the  finite-dimensional  sub- 
h+0  n 
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space 8 

Yh  -  P^Y  C  Y,  Wh  -  (Q|W)_1Yh  C  W,  xh  "  wh  •  Z  C.  X 
are  well-defined  and  we  may  consider  the  approximating  problem 

Vx)  -  y0h*  x  £  R(V»  y0h  -  Vo 

(6.1) 

ph:  shc  \  *  V  ph<x)  ■  php(x)*  x  E  Sh  «  s  n  Xh. 

In  order  to  compare  the  solution  manifolds  of  these  discretized 
problems  with  M(b)  we  need  to  assume  that  the  approximation 
Is  sufficiently  close  to  ensure  that  the  same  local  para- 
metrlzatlon  can  be  used 
on  these  manifolds, 

(see  Figure  7).  As 
usual,  a  stability  con¬ 
dition  Is  required  to 
ensure  the  convergence. 

Without  elaborating 
upon  the  specific  de¬ 
tails  of  this  condition, 
we  sketch  here  only  the 
general  form  of  the  re¬ 
sulting  theorems  (see 
[12]  and  [13]): 

Theorem  6.2;  Under  the  stated  conditions  about  F  and  M(b), 
let  {V,A,T>  be  a  local  parametrization  of  M(b)  at 
xQ  e  M(b)  and  suppose  that  the  Indicated  stability  condition 
holds.  Then,  for  all  sufficiently  small  h  >  0,  the  approxi- 


Flgure  7 


26 


mate  problems  (6.1)  possess  solutions  xQh  e  R(Fh)  for  which 
h-0  oh  0 

Theorem  6.3:  In  addition  to  the  conditions  of  Theorem  6.2, 

let  DF  be  Llpschitz-contlnuous  on  bounded  subsets  of  SCX. 

Denote  by  x:  B  c  T  S  the  representation  of  M  (b)  near  xQ 

defined  on  the  open  ball  B  c.  T,  0  e  B.  Then  there  exists  a 

closed  ball  B  CB.  0  e  B  such  that  for  all  sufficiently 
o  o 

small  h  >  0  the  approximate  problems  have  solutions 

x.  :  B  -*■  S.  for  which 
ho  h 


(6.2) 


| |x(t)-xh(t)| |  <  C| |(I-Ph)Qx(t)| |,  t  e  Bq 

where  the  constant  C  is  independent  of  h  and  t. 

As  a  typical  example  consider  the  mildly  nonlinear 
problem  in  the  weak  form 

a(u,w)  +  |  g(u,X)wd£  -  0,  w  e  H1^),  X  e  Rm, 


where  D  is  a  suitable  bounded  domain  in  R  ,  and  the  bilinear 


form 


:u,w)  -  f  l  a „(€)#-#■ 

J  1,J-1  1J  n3 


dC,  u , w  e  H1(D) 


has  sufficiently  smooth  coefficients  and  is  strongly  elliptic. 

Then  the  operator  K:  L2(fl)  ♦  (SI)  defined  by 

9  °1 

a(Ku,w)  ■  (u, w) Q ,  u  e  IT(n),  w  e  HA(fl),  is  compact  and  the 
problem  is  equivalent  with 


u  +  KG(u,X)  ■  y. 


where  G:  H^(G)  x  Rm  L^(fl)  is  defined  by  the  Cr-map  g.  It 
turns  out  that  this  is  a  Fredholm  map  of  index  m  and  that, 

Oi  — . 

with  Q  as  the  natural  projection  from  H  (Jl)  x  r“  onto 
°1 

H  (n),  the  mentioned  stability  condition  is  satisfies  and 
Theorems  6.2  and  6.3  apply. 

On  the  basis  of  a  priori  estimates  of  the  form  (6.2)  we 
can  now  develop  a  posteriori  error  estimates  along  the  same 
lines  as  for  the  linear  problems.  Some  examples  for  this  are 
presented  in  [7]  and,  as  in  the  linear  case,  the  effectivity 
of  these  estimates  turns  out  to  be  consistently  high.  However, 
the  theory  of  these  a  posteriori  estimates  for  nonlinear 
problems  is  not,  as  yet,  very  well  developed  and,  in  particular, 
such  results  as  Theorem  2.3  are  still  lacking.  In  [7]  it  was 
shown  also  that  a  posteriori  estimates  can  be  constructed 
for  buckling  loads  and  the  location  of  limit  points. 


7.  CONTINUATION  PROCESSES 

For  the  computational  analysis  of  the  solution  manifolds 
considered  here  the  basic  methods  are  general  continuation 
processes  for  the  trace  of  specified  paths  on  the  manifold. 

If  the  problem  formulation  Involves  natural  parameters,  then 
these  paths  are  defined  usually  by  combinations  of  the  para¬ 
meter  values  with  one  degree  of  freedom.  But  interest  may 
center  also  on  paths  which  are  specified  implicitly,  for 
example,  as  paths  in  the  critical  boundary  of  the  problem  (see 
eg.  [17]). 

The  literature  on  continuation  methods  is  rather  large. 
For  a  survey  relating  to  structural  mechanics  we  refer  only  to 


[20]  where  also  further  references  may  be  found.  No  attempt 
shall  be  made  to  present  here  the  details  of  these  methods. 
Broadly  speaking  a  typical  such  process  (PITCON,  [19])  In¬ 
volves  the  following  steps: 

(1)  Compute  the  (normalized)  tangent  vector  at  the 
current  point  on  the  path. 

(2)  Determine  a  new  local  parametrlzatlon  of  the  path 
near  this  point. 

(3)  Determine  a  steplength. 

(4)  Use  the  steplength  to  compute  a  predicted  point 
further  along  the  path. 

(5)  Start  a  corrector  Iteration  from  the  predicted  point 
to  obtain  a  new  point  (approximately)  on  the  path. 

The  overall  process  uses  extensive  feedback  as  the  basis  of 
the  various  decisions.  This  Is  particularly  true  for  the 
above  steps  (2)f  (3),  and  (5).  The  feedback  mechanisms  used 
here  are  often  rather  complex  In  nature,  and  this  raises  the 
question  whether  the  overall  process  can  be  shown  to  be 
adaptive  In  the  sense  of  Section  2.  Up  to  now,  no  such  re¬ 
sults  have  been  formulated,  although,  there  appears  to  be  the 
possibility  of  proving  at  least  some  partial  results.  There 
Is  certainly  a  need  for  more  studies  of  this  area. 

On  the  basis  of  the  estimates  of  the  previous  section  it 
turns  out  that  the  discretization  errors  along  any  continuation 
path  often  vary  considerably.  On  the  other  hand,  our 
objectives  usually  require  that  these  errors  be  maintained 
reasonably  well  within  a  tolerance  Interval.  This  suggests 
again  the  use  of  mesh-reflnements  (and  de-reflnements)  In 


conjunction  with  the  continuation  process. 

In  its  simplest  form  such  a  combined  process  begins  with 
a  mesh  that  meets  a  prescribed  error-tolerance  requirement. 
Then,  during  the  continuation  process,  the  error  estimators 
are  monitored  and,  at  any  computed  point  where  these  esti¬ 
mators  exceed  the  upper  tolerance  a  mesh-refinement  strategy 
is  applied  to  meet  the  tolerance  requirement.  Similarly,  we 
may  work  also  with  a  lower  error  tolerance  and  de-refine  the 
mesh  whenever  the  estimates  fall  below  that  lower  value. 

Figure  8  indicates  the  effect  of  such  a  strategy  for  a 
problem  involving  a  nonlinear  rod.  We  shall  not  present  here 
the  details  of  the  problem-formulation  but  refer,  for  that,  to 
[7].  In  brief,  the  indicated  values  of  N  refer  to  the 
number  of  quadratic  C°-elements  that  were  used  during  the 
corresponding  segments  of  the  path.  The  error  tolerance  was 
2% .  Since  the  continuation  process  takes  relatively  large 
steps  along  the  path,  there  is  no  guarantee  that  between  these 
points  the  error  tolerance  does  not  exceed  the  required 
tolerance.  A  more  sophisticated  control-strategy  may 
alleviate  the  situation  somewhat.  At  the  same  time,  problems 
of  this  type  cause  difficulties  in  developing  theoretical  re¬ 
sults  of  the  type  of  Theorems  2.1  and  2.2  for  the  combined 
continuation  and  mesh-refinement  processes.  The  level  of  ex¬ 
pected  difficulty  of  such  results  for  nonlinear  problems 
suggests  that  at  first  some  emphasis  should  be  placed  upon 
carefully  designed  experiments  which  confirm  or  refute  various 
conjectures  about  the  adaptivity  of  the  feedback  mechanisms 
that  are  used  in  these  settings. 


Figure  8 
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